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Abstract. 

Properties of hcp-Ti such as elastic constants, stacking faults and 7— surfaces 
are computed using Density Functional Theory (DFT) and two central force 
Embedded Atom interaction Models (EAM) [T] [2]- The results are compared 
to previously published calculations and to predicting models. Their implications 
on the plastic properties of hcp-Ti are discussed. 



PACS numbers: 62.20.-x Mechanical properties of solids, 71. 15. Mb Density 
functional theory, local density approximation, gradient and other corrections, 
61.72.Lk Linear defects: dislocations, disclinations 
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1. Introduction 

The plastic behavior of hexagonal compact metals (about twenty) is dominated by the 
movement of dislocations with the shortest Burgers vector: #(1120) (see D. Caillard 
and J.-L. Martin [3J, and M. H. Yoo and coworkers [I] for recent reviews). For example, 
most of the transition metals have a basal slip plane, but three of them (Zr, Ti and Hf ) 
and a number of rare-earth metals (Gd, Ru, Tb, Hf, Dy and Er) present a dominant 
prismatic slip plane 4, 5, 6, 7\. The intrinsic characteristics of dislocations (structure, 
stability, mobility, formation, multiplication) as well as their extrinsic characteristics 
(interaction with impurities, point defects, grain boundaries...) are therefore essential 
ingredients for a good theoretical description of the plastic behaviour of hexagonal 
close packed (hep) metals. 

The core structure of the #(1120) screw dislocation in hep metals has been 
subjected to only a few theoretical studies (see D. J. Bacon and V. Vitek for a 
recent review [B]). B. Legrand, who introduced an explicit treatment of the electronic 
degrees of freedom in a tight-binding (TB) approach [7J [5], was the first one to 
obtain a prismatic spreading in an hep metal (Ti). Moreover, this prismatic spreading 
was energetically more favorable than the basal one by 24 meV/A. He argued that 
this prismatic preference was related to an increased basal stacking fault energy of 
transition metals with d fillings between 1.5 and 2.5 electrons. This effect results 
from the directional <i-covalent bonding of partially filled d bands, which cannot 
be described with pair potentials nor with more sophisticated central force N-body 
potentials. More precisely, he established [3 [10] on the basis of his TB calculations 
(and pseudopotential calculations for divalent and trivalent metals), a clear correlation 
between the basal or prismatic slip in hep metals and the ratio R = (713/044)/ {l-p/G^o)- 
7b and 7 p are respectively the basal I2 and the prismatic stacking faults excess energies 
and C44 and C66 the shear elastic constants governing the shear deformation leading to 
the corresponding stacking faults. This ratio measures the relative facility to form both 
stacking faults and thus the easiness of dislocation dissociation in the corresponding 
planes. For R > 1, the observed slip should be prismatic whereas it should be basal 
for R < 1. One must quote that, in all his calculations, he never found a stable 
prismatic stacking fault, but the 7 P -surface presented a deep valley in the (1210) 
direction. He defined then the stacking fault excess energy as an average value of 
the 7 P surface energy in this direction, corresponding to a continuous dislocation 
spreading. Notwithstanding, in principle, their unability to describe the directional 
character of the bonding in Ti, pair potentials and N-body central force potentials 
have been able to reproduce the preferential prismatic core structure of the |(1120) 
screw dislocation in some cases [IT] [T2] . 

More recently, A. Girshick and coworkers [133 derived a bond order (BO) potential 
for Ti [14] allowing them to deduce more accurate values of the excess energies of the 
basal and prismatic dislocation cores. With the BO potential, in agreement with B. 
Legrand [TJ [8] , they found that the prismatic core was the most stable one (by about 20 
meV/A relative to the basal one). At odds with the BO potential, the N-body Finnis- 
Sinclair central force potential |15j that they used for comparison, favored the basal 
spreading (by about 25 meV/A). In agreement with B. Legrand, they did not found any 
stable prismatic stacking fault. Most importantly, from their BO calculations, they 
obtained a R = 0.4 value, which means a preferential basal spreading in contradiction 
with their determination of the core structure energies, which questions the validity 
of Legrand's criterion. 
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Up to now, DFT calculations on the \ (1120) screw dislocation in hep metals 
(Zr and Ti) were reported by C. Domain and A. Legris [TBI E] on very small 
clusters and by N. Tarrat et al. [IB] on slightly larger ones. Both studies showed 
a prismatic spreading with a core structure in overall agreement with the previous TB 
calculations [8j [13] . Very recently, DFT calculations on very large clusters using 
flexible boundary conditions have been reported [19] . They found similar prismatic 
spreadings for the ^(1120) screw dislocation than the ones obtained by Tarrat et al., 
but argued that the prismatic core structures are slightly different and that one of 
them is actually a metastable one. 

The purpose of this paper is therefore to present accurate calculations of the 
bulk properties, stacking faults and 7-surfaces which are essential properties to gain 
a deeper understanding of the relationship between the electronic structure of hep 
Ti and its plastic behavior. The DFT results are compared to previously published 
calculations (TB 0[H], BO potential [T3]. MEAM potential [3D] and DFT calculations 
[ini [TT] [5UJ ) and to similar calculations performed with recently developed accurate 
semi-empirical EAM potentials [TJ 2 . These EAM potentials were specifically designed 
to reproduce a large data base of properties in TiAl [JJ or Ti [2] , either experimentally 
measured or deduced from ab initio DFT calculations. These data bases include a 
large panel of configurations far from the ideal hep lattice, so these potentials might 
be accurate for the description of the dislocation core structure. 

After giving the simulation details in Sec. [2] we present the results obtained on 
the bulk properties, stacking faults and 7-surfaces of hep Ti using two EAM potentials 
and DFT calculations in Sees. 13. 1[ 13. 21 and [331 The obtained results are discussed and 
compared to the literature. In particular their implication on the plastic behavior of 
hep Ti is examined. 

2. Simulation details 

2.1. DFT calculations 

The DFT calculations have been performed with the SIESTA |2T] code in which the 
orbitals are developed on a local basis set. The PBE-GGA gradient approximation [22] 
has been used for the exchange and correlation functional and the pseudopotential was 
a Trouiller-Martins [33J one, with the 3p4s3d states as valence states. A polarized 
doublc-C basis set was employed for the 4s electrons and a single-C basis set was 
employed for the 3p and 4d electrons. 

A 500 Ry real space grid cutoff and a Methfessel-Paxton smearing of order one 
with an electronic temperature of 300 meV were used in all calculations. For the bulk 
calculations, performed with the conventional primitive cell, a 11x11x7 k-point mesh 
was used for the first Brillouin zone sampling, which was extended to a 16x16x12 
mesh for the determination of some of the elastic constants. For the other calculations, 
the specific k-point-mesh will be given in the corresponding section. 

2.2. Semiempirical potential calculations 

The two semiempirical potentials were of the embedded atom method [24j type. The 
first one (referred here after as ZM) was developed by R. R. Zope and Y. Mishin [JJ for 
the TiAl alloy. They fitted their Ti potential to experimental data (lattice parameters: 
a, c/a, cohesive energy E c and the five independant elastic constants C^) and to the 
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ab initio DFT volume-pressure curves of various crystal structures (hep, face-centered 
cubic (fee), body-centered cubic (bec), simple cubic (sc) and ui). Its cutoff distance is 
in between the fourth and fifth neighbour shells. The second one (referred here after 
as HKV) was developed by Hammerschmidt and coworkers [5] for the description of 
grain boundaries. It was also fitted to bulk Ti properties and to a large data base 
of low coordinated configurations, namely small clusters (up to 8 atoms), surfaces 
and adatom surface diffusion. Its cutoff distance is in between the third and fourth 
neighbour shells. Both of them are then expected to be able to describe configurations 
far from the ideal hep structure. 

In all calculations where periodic conditions were applied in a given direction, the 
simulation cell dimension in that direction was set to at least five times the interaction 
cutoff in order to avoid any spurious interaction of an atom with its periodic images. 

3. Results 

3.1. Bulk properties 

In Tables [1] and [H the bulk properties of hep Ti calculated with DFT and the two 
EAM potentials are compared to previous calculations performed with DFT the 
BO potential [TJ] and a new MEAM potential [3D] and to experimental data. 

Table 1. Properties of bulk hep Ti predicted by the DFT calculations and the 
two EAM potentials [TJ [2] compared to previously published calculations 1 1 41 1 20 1 
and to experimental data I25| . Lattice parameter a is in A, Vo is in A 3 /at., the 
bulk modulus B is in GPa. 





a 


c/a 


Vo 


B 


DFT (this work) 


2.996 


1.588 


18.49 


110.2 


DFT [15] 


2.947 


1.583 


17.55 




ZM (this work) 


2.951 


1.585 


17.64 


110.5 


HKV (this work) 


2.969 


1.590 


18.02 


110.4 


BO[H] 


2.950 


1.587 


17.65 


113.6 


MEAM [20] 


2.930 


1.596 


17.40 




Exp. [25, 


2.951 


1.588 


17.65 


110 



The equilibrium volume and bulk modulus are obtained with a similar accuracy 
in all cases. For the elastic constants, the agreement between the DFT calculations 
and experiment is very satisfactory (deviations less than 10%), better than with the 
EAM, BO or MEAM potentials (deviations up to ~ 40%, 20% and 15 % for the 
EAM, the BO and the MEAM potentials, respectively). Notably, there is a large 
deviation on Cg6 with the EAM potentials. It is particularly interesting to note 
that the Cauchy pressures, which are due to N-body interactions only, are rather 
well reproduced by the present DFT calculations and the MEAM ones, whereas the 
deviations from experimental values are large (up to 90%) using the EAM or the BO 
potentials for which C13 - C44 is overestimated, whereas C12 - C66 is underestimated. 
The MEAM potential appears then as the best interaction model for the description of 
the Ti elastic properties. The small differences between the present DFT calculations 
and the ones of Ref.[20] are certainly due to the use of a different GGA functional and 
a different basis set. 
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Table 2. Elastic constants of bulk hep Ti predicted by the DFT calculations 
and the two EAM potentials [T][2] compared to previously published calculations 
1 141 1201 and to experimental data |25| . All elastic constants Gy and Cauchy 
pressures (CPi=Ci2-C66 and CP2=Ci3-C44) are in GPa. Note that the bond 
order potential of Ref. |14| was fitted exactly on the Cn, C33 and C44 elastic 
constants. 





C n 


C12 


C13 


C33 


C44 


^66 


CPi 


CP 2 


DFT (this work) 


183.4 


84.6 


63.8 


204.9 


48.8 


49.4 


35.2 


15.0 


DFT|20] 


172 


82 


75 


190 


45 


45 


37 


30 


ZM (this work) 


186.2 


69.5 


76.2 


189.4 


46.4 


58.3 


11.2 


29.8 


HKV (this work) 


188.6 


65.4 


67.4 


216.9 


45.8 


61.6 


3.8 


21.6 


BO[H] 


176.0 


74.0 


83.3 


190.5 


50.8 


51.1 


22.9 


32.5 


ME AM [20] 


174 


95 


72 


188 


58 


39.5 


55.5 


14 


Exp. 25 


176.1 


86.9 


68.3 


190.5 


50.8 


44.6 


42.3 


17.5 



3.2. Stacking faults 

The excess energies of the two intrinsic Ii (ABAB|CBCB) and I 2 (ABAB|CACA), 
and the extrinsic Ie (ABAB|C|ABAB) elementary stacking faults in the basal plane 
have been determined, as well as the easy prismatic stacking fault. They have been 
calculated in a slab geometry with free surfaces. 
The excess fault energy 7 is then given by: 

7 = -E-faultcd — -^perfect, (1) 

where -Bf au itcd is the total energy of the slab containing the fault, and -E pe rfcct is 
the total energy of the perfect slab. For the DFT calculations, the size of the vacuum 
(4 interatomic distances, i.e. ss 9.5 A) and the number of planes (14 atomic planes 
for Ii and I 2 and 15 atomic planes for 1^) have been adjusted to have good converged 
results. The slabs were then relaxed, using a mesh of 12x12x2 k-points. In the EAM 
simulations, the slabs were made of 40 atomic planes and were replicated 8 times in 
the fault plane, so that the dimensions of the supercell are about 5 times the cutoff 
radius of the potential. 

The stacking fault energies, together with the (0001) surface energy, are presented 
in Tableland compared to those obtained with the EAM potentials, to experiments 
and to previously published results (7] [8j US] QU OH [H [19] . All the DFT calculations 
give very similar results expect those of Refs.[16l [17] for the prismatic stacking fault 
value which they found 50% lower. This difference comes from the strong dependence 
of the excess energy with the number of planes in the slab. Indeed, in Fig. 13. 2\ the 
evolution of the prismatic stacking fault excess energy as a function of the number 
of atomic planes included in the calculation is presented. The excess energy starts 
to converge when the number of planes is at least greater than 16, but with still 
noticeable oscillations with the number of planes. These oscillations are known to be 
due to possible quantum size effects for specific slab geometries which lead to long 
range interlayer oscillating relaxations [35] [33J [3H [33] [33] and consequently make the 
DFT calculations very sensitive to the number of planes. The present DFT value 
reported in Tab. [3J is the mean value for calculations with a number of planes larger 
than 16. Considering the very indirect experimental determination of the stacking 
fault excess energies, the DFT results are satisfactory. 
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Table 3. Hep Ti (0001) surface energy 7 S and stacking fault excess energies 
(mj/m 2 ): 7i 2 , 7i E (0001) stacking faults and 7 P prismatic stacking fault. 
Present DFT and EAM (ZM T and HKV [2]) potentials calculations are compared 
to experimental data and previous calculations with MEAM, BO, TB and DFT 
methods, 





7s 


Tii 


7l2 


7i E 


Tp-casy 


DFT (this work) 


2048 


148.6 


259.1 


353.1 


250.0 


ZM (this work) 


1263 


30.6 


54.3 


82.4 


364.6 


ZMpQ 


1725 


31 


56 


82 




HKV (this work) 


1188 


33.6 


64.5 


94.2 


353.6 


HKV0 


1185 










A-FSH3I 




33 


64 


94 


253 


Bopg 






110 




260 


TB0HIS] 






290-370 




110 -140 


TB[26] 




44 


118 






DFT [16] H3 






291 




174 


DFT [27] 


2045 




287 






DFT 20 , 19 


1938 




272-320 




220 


MEAM [ID] 


1474 




170-172 




297 


Exp. 


2100128 , 1920[29 




> 300 30 




150[31 



The two EAM potentials give identical results but those results are quite different 
from the DFT values in spite of the fact that they were chosen because they were 
developed in order to reproduce defects with high angular variations in titanium. 
These differences may have important consequences for dislocation core structure 
calculations as well as for the evaluation of the Peierls stress. The results obtained 
in the TB scheme or with the BO or MEAM potentials are much closer to the DFT 
ones. This can be explained by a better reproduction of the N-body interactions. 
However, surprisingly, some of those TB or BO results [HE] [8] reproduce properly the 
basal stacking fault, and noorlv the nrismatic. one. whereas some others |14U26j do the 
reverse. 

300 

280 
260 

^ 240 

200 
180 
160 

4 8 12 16 20 24 28 32 36 40 44 48 52 56 60 

Number of atomic planes 

Figure 1. Evolution of the DFT prismatic stacking fault excess energy as a 
function of the number of atomic planes. 
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The MEAM potential still appears as the best approximated interaction model, 
however it is still unable to reproduce the correct energetic ordering between the basal 
I2 and the easy prismatic stacking faults which could lead to an erroneous preferential 
basal spreading of the |(1120) screw dislocation and consequently a poor description 
of the plastic properties. 



3.3. ^-surfaces 

For the determination of the 7— surfaces, similar systems have been used, i.e. slabs 
with free surfaces, of 14 atomic planes for the basal surface and of 16 atomic planes 
for the prismatic one. In the EAM simulations, the slabs were made of 40 planes of 
64 atoms in the basal case and of 40 planes of 70 atoms in the prismatic case, leading 
to systems of 2560 and 2800 atoms, respectively. 

The slabs were cut into two parts, which were then shifted with respect to each 
other by steps of 0.005 x a in the empirical simulations and of 0.01 x a in the DFT 
ones. The atoms were allowed to relax only in the direction perpendicular to the slab 
surface. 

The basal 7— surfaces obtained with the two EAM potentials have a similar global 
shape than the one obtained in DFT (Fig. [2j left). However one can note differences 
in particular on the basins, which are quite deep using empirical potentials, and are 
almost inexistant using DFT. These basins correspond to the basal stacking faults 
(~fI E in Table [3]) , which were indeed found much lower than the DFT one. The only 
difference between the two EAM potentials appears to be a more flat maximum in the 
HKV potential case than in the ZM case. 

The two EAM prismatic surfaces (Fig. [5J right) present only small differences. 
The HKV surface is slightly flatter than the ZM one. However their differences with 
the DFT one are more important. The maxima are quite lower and the basins are 
deeper. One can also see that the DFT 7— line along a | (1120) direction differs notably 
from those of the two EAM potentials. This result is important since it is commonly 
accepted that the shape of the 7— line is correlated with the dislocation gliding. So 
an overestimation of the excess energy can lead to an error in the spreading plane. 
This can be more easily seen in Fig. [3] where the basal and prismatic 7— lines along 
a |(1120) Burgers vector are depicted. In the basal case (Fig. [3Ka)), the differences 
between the DFT line and the EAM potentials ones are negligeable. This is not true in 
the prismatic case (Fig. [3Jb)) where the DFT line exhibits a local minimum contrary 
to the two EAM potentials and correspondingly its maximum is much lower. In the 
EAM calculations, there is also a local minimum but it is slightly shifted from the 
|(1120) line (Fig. [2j right). This local minimum of ~ 18.6 mJ/m 2 at the center 
of the line in the DFT case indicates the existence of a stable stacking fault in the 
prismatic plane, and thus a possible dissociation into partial dislocations. This result 
was already obtained in Zr QH [T7] and Ti QH [TTJ HI] by DFT calculations but not 
with empirical or semi-empirical potentials [131 119] . 

4. Discussion and conclusion 

A detailed study of bulk properties of hcp-Ti using the state of the art DFT scheme or 
approximated central force EAM interaction models has been presented. Important 
properties for the description of dislocations have been studied: the elastic constants, 
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HKV potential 




ZM potential 




Figure 2. Hep Ti gamma surfaces. Left: basal plane, right: prismatic plane. Top 
and middle: HKV and ZM EAM potentials calculations, respectively. Bottom: 
DFT calculations. 




Prismatic line 
400 1 




a/3 < 1 1 -2 > 



Figure 3. Hep Ti basal (a) and prismatic (b) 7— lines in the ^(1120) Burgers 
vector direction. Circles: DFT calculations. Squares and rhombus: ZM and HKV 
EAM potentials calculations, respectively. 



the stacking faults and the 7— surfaces. Using the calculated stacking fault energies 
one can evaluate the R ratio of B. Legrand [5] for the determination of the easy slip 
plane for the |(1120) screw dislocation (see the Introduction). 

The DFT results give R = 1.1, a value slightly higher than 1, leading to a 
preferential prismatic spreading as experimentally observed, and in agreement with 
our previous DFT calculations of the 1(1120) dislocation core structure [18]. Our 
DFT R value is consistent with the one deduced from previously published DFT 
calculations (JTH1 El [201 H]) which always lead to R values slightly larger than 1. 
The DFT ratio, just above 1, is much less than the one initially obtained by Legrand 
(2.5) in his TB calculations. His higher value was mainly due to the underestimation 
of the prismatic 7 Pi easy stacking fault energy [8] by a factor of about 2. Unexpectedly, 
in its TB - BO calculations, Girschich and coworkers [13] found R = 0.4, a value much 
lower than 1, which is supposed to lead to a preferential basal spreading, whereas they 
observed a prismatic one. In that case, surprisingly the difference comes essentially 
from an understimation of the basal stacking fault energy by a factor of about 2 [13] . 
A similar underestimation of this basal stacking fault was also obtained by Bere [26] . 
It is surprising that these three TB calculations give so different values for the stacking 
fault energies, in particular for the basal one. Indeed in the TB approach, to obtain 
the right energetic ordering of the hep, fee and bec phases in compact close-packed 
structures, which conditions a good value for the basal stacking fault, one needs to 
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make a calculation involving up to the sixth moment of the electronic density of states 
for the electronic energy calculation with the used recursion method. In Legrand's 
calculations six moments were used, whereas in Girshick's and Bere's calculations 
9 and 32 moments were used respectively. So, in these TB calculations the basal 
stacking fault excess energy, which originates mainly from the electronic part of the 
total energyJJ should be described with increasing accuracy with the number of used 
moments. As expected the basal stacking fault energy is rather well evaluated by 
Legrand, but unexpectedly this evaluation is bad in Girshick's and Bere's calculations 
with much too low values. So, increasing the accuracy in the description of the 
electronic band structure, instead of improving the description of the basal stacking 
fault, makes it dramatically worth. At odds, it seems to improve the description of 
the prismatic stacking fault excess energy. In that case, however, since the interatomic 
distances are quite modified, the excess energy will not be due to the sole electronic 
energy. 

Using the EAM results, we obtain a 0.2 R value, leading to a preferential basal 
spreading as we do have observed in our EAM calculations of the 1(1120) dislocation 
core structure [T5] . The introduction of angular forces in the MEAM model is supposed 
to allow a better description of the basal stacking fault. Indeed, the calculated MEAM 
[2"UI [T9"] basal ji 2 excess energy is three times higher than the EAM one. However, it 
is still lower than the DFT one leading to a R ratio significantly lower than 1 (.4), 
since the prismatic fault energy is properly evaluated. 

As previously discussed by Domain [17] . the Legrand criterion is too simplified. 
In particular, it does not take into account the energy cost associated to the edge 
components of the dissociated basal dislocation. This extra energy cost means a 
preferential prismatic spreading even with R values lower than one. This explains 
the preferential prismatic spreading observed in the DFT calculations even with an 
R ratio only slightly larger than one, and similarly could explain the BO observed 
preferential spreading with R = 0.4. More puzzeling is the large observed dispersion 
in the TB based calculations of the basal stacking fault and its still poor description 
with the MEAM potential, since it may lead to erroneous dislocations calculations in 
these approaches. These results suggest that only a DFT approach or a more precise 
TB approach, like the DFT-based tight-binding approach (DFTB) should be able to 
properly describe the dislocation core structures in hep Ti and consequently its plastic 
behavior. 
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| The interatomic distances and the coordination numbers are unchanged up to the second neighbor 
shell (included), which means an essentially unchanged repulsive energy in the model. 
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